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We examine tlie ratio of tlie decay rate of tlie r/c into light hadrons to the 
decay rate into photons and find that most of the large next-to-leading-order 
(NLO) correction is associated with running of the strong coupling a^. We re- 
sum such contributions by analyzing final-state chains of vacuum-polarization 
bubbles. We show that the nonperturbative parts of the bubble chains can 
be absorbed into a color-octet matrix element, once one has used contour 
deformations of the phase-space integrals to cancel certain contributions. We 
argue that these contributions are incompatible with the uncertainty princi- 
ple. We also argue that perturbation theory is reliable only if one carries out 
the phase-space integrations before the perturbation summation. Our results 
are in good agreement with experiment and differ considerably from those 



that one obtains by applying the scale-setting method of Brodsky, Lepage, 

and Mackenzie to the NLO result. 
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I. INTRODUCTION 



In this paper we consider the decay of the heavy quarkonium 7]^ into hght hadrons. 
According to the analysis of Ref. the decay rate is given by the Nonrelativistic Quantum 
Chromodynamics (NRQCD) factorization formula 

2 Im f 

nvc-^m = E^{vc\oM. (1.1) 

Here, rric is the charm-quark mass, the On are operators in NRQCD with mass dimension 
dn, and 2Im/„/m^"~^ is a short-distance coefficient. A similar formula applies for the 
electromagnetic decay rjc 77. 

In general, there are difficulties in obtaining accurate results from the factorization for- 
mula (|1 ■ 1| ) ■ The sum is an expansion in powers of the heavy-quark-antiquark relative velocity 
V. Only the first few operator matrix elements in that sum are known at all, and they are 
not known very accurately. The short- distance coefficients may, in principle, be computed 
in perturbation theory. However, the perturbation series is only asymptotic. In particular, 
it exhibits a factorial growth that is associated with the presence of renormalon singularities 
in its Borel transform P|. Such singularities give rise to ambiguities in the evaluation of the 
perturbation series. 

One can avoid some of these difficulties by considering the ratio of the rate for rjc into 
light hadrons to the rate for rjc into two photons. We denote this ratio by R. In R, the 
dependences of the decay rates on the operator matrix elements cancel in leading order 
in V, and the corrections of relative order also cancel These cancellations occur 
because, aside from an overall color factor, the tree-level Feynman diagrams for QQ —>■ gg 
and QQ '-ff are the same. One consequence of these cancellations is that renormalon 
ambiguities associated with initial-state virtual-gluon corrections also cancel at the 
level of accuracy of relative order v"^. However, there is still the potential for renormalon 
ambiguities associated with final-state corrections to appear. 

Through next-to- leading order in f^, the ratio R is given by the ratio of the short- distance 
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coefficients of the leading-order operators Through next-to-leading order (NLO) in as, 
one obtains 



where 



' /199 137r2 8 \ aJ/i) ^ , , /i^ 



:i.2) 



mc is the pole mass of the charm quark, is the number of light quarks (nj = 3 for 
charmonium), [3q is the leading coefficient of the beta function, with [3q = (33 — 2nj)/(127r), 
and, in QCD, A^, = 3, and Cp = 4/3. 

Setting the factorization scale fj, equal to 2mc and taking as{2mc) = 0.247, one finds that 

i?^^°(2me) = 2.1 X lOl (1.4) 

This is to be compared with the world-average experimental value 0] 

^Exp ^ 3 3 ^ ^ 3 ^ ^g3 _ ^^ 5^) 

However, the next-to-leading-order QCD correction is quite large. For as{2mc) = 0.247, it 
is about 1.1 times the leading-order term. Hence, there is some question as to the reliability 
of the perturbation expansion. 

Furthermore, the dependence on the factorization scale is large. For example, taking 
/i = rric, with as{mc) = 0.350, we obtain 

i?NLO(^j ^ 5 ^ ^03^ 

which is factor of 2.4 larger than R^^'^ {2mc) . 

One popular choice of scale is that given by the Brodsky-Lepage-Mackenzie (BLM) 
method In this approach, all terms proportional to Uf are absorbed into the running 
coupling as{fi). For R, one obtains /iBLM ~ 0.52mc This yields a next-to- leading-order 
coefficient that is only about 0.5 times the leading-order coefficient. However, 



i^'^'^^l/iBLM) = 9.9 X 10^ (1.7) 

which is much larger than either the theoretical result at /i = 2mc or the experimental value. 

We can gain some insight into the origins of the large next-to-leading-order correction 
and the strong scale dependence by re-writing Eq. ( |1 . 2|) as follows: 



i?^L°(/i) = /?o(/x)<il + 



'37 13n^\ „ (IQ _ /i^ 



+ 7r/9o — + 21n 



TT 



1.8) 



2 8 / ' \?, Aml^ 

In the next-to-leading-order correction, the first term in square brackets has a value of about 
2.46, while the term 167r/5o/3 has a value of 12. Hence, we see that most of the large next-to- 
leading-order correction is proportional to j3o. That is, it arises from the one- loop corrections 
associated with the running of a^. 

This suggests that we control the large next-to-leading-order correction by resumming, to 
all orders in perturbation theory, the contributions associated with the one-loop running of 
Us- Such a resummation would have the added benefit of reducing the scale sensitivity, since 
the resummed contribution would be invariant under a change of scale, aside from two-loop 
corrections to the running of a^- Unfortunately, there is some arbitrariness in the resum- 
mation procedure. The contributions to the running of that arise from quark loops are 
contained entirely in chains consisting of quark vacuum-polarization bubbles connected by 
gluon propagators. However, the contributions that arise from gluon loops may, depending 
on the gauge, be associated with diagrams other the gluon vacuum polarization. 

One unambiguous procedure for carrying out such a resummation is the method of 
"naive non-Abelianization" (NNA) |]^ . In this method, one resums the chains of quark-loop 
vacuum-polarization bubbles and then accounts for the gluonic contributions by promoting 
the quark-loop contribution to (3q to the full QCD (3q. This method is a generalization of the 
BLM scale-setting method. However, as we shall see, it leads to a very different numerical 
result. An alternative method for carrying out the resummation is to resum both the quark- 
loop and gluon-loop vacuum-polarization bubble chains in the background-field gauge [jlOl . 



In this gauge, all of the corrections associated with the running of as are contained in the 
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quark and gluon vacuum-polarization diagrams. It turns out that the background-field-gauge 
method gives numerical results that are close to those of the NNA method. 

Vacuum-polarization insertions first appear, in relative order as, in the final-state gluon 
legs. Therefore, in order to resum the leading bubble-chain corrections, we need consider only 
vacuum-polarization bubble chains in the two final-state gluon legs. By "leading bubble- 
chain corrections," we mean those corrections that contain one power of the large coefficient 
in the second term in square brackets in Eq. (|1.8| ) for each power of as- 

The integrals over the final-state phase space include regions in which the constituents 
of the vacuum-polarization diagrams (quarks or gluons) have a small invariant mass. Hence, 
the calculation is, in principle, sensitive to nonperturbative contributions. In the bubble- 
chain series, this sensitivity may manifest itself as a factorial growth that is characteristic 
of the presence of renormalon singularities in the Borel transform of the series. 

In the case of a single chain of vacuum-polarization bubbles, we show, by making use of a 
contour-deformation argument in the complex invariant- mass plane, that the low-virtuality 
contributions, and, hence, the factorial growth, cancel when one carries out the phase- 
space integration before carrying out the perturbation summation. Somewhat surprisingly, 
even when factorial growth in the series cancels, there is an ambiguity that arises in the 
bubble-chain series: one obtains different results, depending on whether one carries out the 
perturbation summation before the phase-space integration, or vice versa. We argue that 
the perturbation series is reliable only when one carries out the phase-space integration first. 

For the more general case of two bubble chains, there is a true sensitivity to the low- 
virtuality region, even after contour-deformation arguments have been applied. Indeed, the 
perturbation series exhibits a factorial growth that is characteristic of the presence of renor- 
malon singularities in its Borel transform. We show that the low-virtuality contributions 
can be identified with contributions to the decay rate that arise from the matrix element of 
a color-octet NRQCD operator. Hence, by taking into account the color-octet contribution, 
one can systematically remove the low-virtuality region from the perturbative calculation. 
In this analysis, contour deformation arguments are crucial to demonstrate the cancellation 
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of certain low-virtuality contributions that do not correspond to matrix elements of NRQCD 
operators. We argue that such contributions are inconsistent with the uncertainty principle. 

The color-octet matrix element is inherently nonperturbative in nature. We do not 
calculate its contribution, which is of relative order v^. Instead, we estimate the size of the 
contribution and treat it as an uncertainty in the calculation. 

The remainder of this paper is organized as follows. In Sec. ||, we compute the con- 
tributions of the vacuum-polarization bubble chains to all orders in a.^- We re- arrange the 



real and virtual contributions to obtain manifestly infrared finite expressions. In Sec. |T|, 
we discuss the contributions from the nonperturbative regions of small invariant mass, first 
for a single bubble chain and then for two bubble chains. We identify the nonperturbative 
contributions that remain after the contour-deformation cancellation with contributions to 



the matrix elements of a color-octet operator. Sec. |^ contains a re-arrangement of the 
bubble-chain contributions that is efficient for numerical evaluation. In Sec. 0, we present 
our numerical results, and in Sec. |V| we give a summary and a discussion of our findings. 

II. THE BUBBLE-CHAIN CONTRIBUTIONS TO R 

In this section, we compute the contributions to R from the insertion of any number 
of vacuum-polarization bubbles into the two final-state gluons in the lowest-order Feynman 
diagrams for QQ gg. 

These contributions are obtained by cutting the diagrams for the heavy-quark-antiquark 
forward scattering amplitude through light-quark and gluon final states in all possible ways. 
A cut can pass through two vacuum-polarization bubbles, a final-state gluon and a vacuum- 
polarization bubble, or two final-state gluons. We call such contributions real-real, virtual- 
real, and virtual-virtual, respectively. These contributions are separately infrared divergent. 



but the Kinoshita-Lee-Nauenberg (KLN) theorem [^] ensures that the infrared divergences 
cancel in the complete expression for the decay width. 
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A. The Effects of Bubble Chains in the Hadronic Cross Section 



In NRQCD, at leading order in f, the decay of a QQ ^Sq state proceeds through the 
color-singlet operator 

0^{}S^)=^\xH, (2.1) 

where ip is the two-component Pauli spinor that annihilates the charm quark and is 
the two-component Pauli spinor that annihilates the charm antiquark. According to the 
NRQCD factorization formula the contribution of this operator to the decay rate in 
NRQCD is given by 

7i('^o) 



T{H) = 2Im 



ml 



{H\OiCSo)\H), (2.2) 



where if is a hadronic state. 

We can compute the short-distance coefficient 2 Im [/i(^S'o)/m^] in perturbation theory 
by taking the state if to be a free QQ state. Then, Eq. (|2.2| ) becomes T{QQ), the rate 
for the annihilation of a free heavy-quark-antiquark pair into gluons through the operator 

At leading order in the heavy-quark-antiquark velocity v in the quark-antiquark center- 
of-mass frame, the amplitude for the decay of a heavy quark and antiquark in a color-singlet 
^5*0 state with momenta pi and p2, respectively, into two massive gluons with momenta k 
and /, polarizations fi and z/, and color indices a and b, respectively, is given by 

^S(QQ) = ^TV(T.T.) 2V2g' ^^^^ . (2.3) 

where g is the strong coupling constant, Ta is an SU{3) color matrix in the fundamental 
representation, and the trace is over the indices of the SU (3) matrices. In deriving Eq. ( p. 3D , 
we have made use of the spin-0 projector, accurate to leading order in v. 



J2 u{pi,Si)v{-p2,S2){s = 0\SiS2) = —y= 



1 / y , N 1 + 70 



2 



i + mj — - — 75(^2 -"^c)- (2.4) 
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Taking the absolute square of the quantity in Eq. ( p.3|) , inserting vacuum-polarization 
bubbles in the gluon propagators, dividing by two for Bose statistics, and integrating over 
the phase space, we obtain the bubble-chain contribution to the QQ decay width: 



x2Im[K(™)(a;)]2Im[ir(")(y)](27r)^5^(pi +P2-k-l), (2.5) 



where 



x= k'^/{Aml), (2.6a) 
y= t^/iAml). (2.6b) 

Aside from a polarization-tensor factor, 

iK^^\x) = Hn(x)r (2.7a) 

4m^(2; + le) 

is the amplitude for the mth-order contribution to the bubble-chain gluon propagator. We 
also define the complete propagator 



K{x) = K^'^\x). (2.7b) 

m=0 

n(x) is defined by 

n^.(fc) = {k^g^u - k^k,)U{x), (2.8) 

where 11^,^ (fc) is the amplitude for a vacuum-polarization bubble with external momentum 
k and polarization indices and u. We note the useful kinematic relations 

ko = mc{l + X — y), (2.9a) 

lo = m,{l-x + y), (2.9b) 

k2 = P = mlil + x'^ + y'^-2x-2y- 2xy). (2.9c) 
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Using Eq. (|2l5| ) and 



(QQ|Oi(^^o)|QQ) = 2iV„ 



(2.10) 



we compare the left and right sides of Eq. (|2.2|) to obtain the short- distance coefficient 



2Im 



^Bub 



2N, 



(2.11) 



Then, taking the hadronic state H in Eq. ( |2.2| ) to be the rjc, we obtain the decay width of 
the rjc into two gluons in the bubble-chain approximation: 



E 



c ^,„=o^ (2vr)4"^-^°^y (27r)4"--^ '^^'^ V2mfco + P, 
X21m[ir(")(x)] 2Im[ir(")(y)](27r)45^(j9i + J92 - A; - /)(r7e|Oi('5o)hc). (2.12) 



dH 



1 



In the rjc width into two gluons ( p.l2|) , we change variables by making use of Eq. ( p. 61) , 
which implies that 2kodko9{ko) = Amldx 9{x) and 2lodlo0{lo) = 4:mldy9{y), and Eq. ( p. 91) . 
We use the energy-momentum-conserving 5-functions to carry out the integrations over the 
3-momenta. Then, dividing by the width into two photons, 

327ra2 



(r/,|Oi(i5o)|r/e), 



(2.13) 



we obtain the contribution to R: 

"1 dx dy 

m,n=0 



i?^"^ = i?o E / /2lm[4mlK^^\x)]2lm[AmlK(-\y)]f{x,y)e{l-V^-^), 
J~7-n Jo 2tx Jo 2tx 

(2.14) 



where 



f{x,y) 



[1 - 2{x + y) + {x - yff/^ 



(2.15) 



{1-x-yy 

In order to compute R^^^ from Eq. (|2.14|) , we need the value of a single vacuum- 
polarization bubble. A standard calculation of the light-quark bubble in dimensional regu- 
larization yields 



U(x) 



-asPo (bx~ 
e ^ 
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(2.16) 



where 



1 + e -7 + In 47r + 



In 



1 + e (^-7 + ln47r + ^ + + 0(e' 



(2.17a) 



(2.17b) 



In Eq. ( p.l6| ), we have renormahzed the ultraviolet divergence by subtracting the ultraviolet 
pole in e and some associated constants, e now plays the role of an infrared regulator (e < 0). 



C is a subtraction-scheme-dependent constant; in MS, C = —5/3. We note that the phase 
of b in Eq. ( p.l7a| ) and the reality of a are exact to all orders in e. For the light-quark 
bubble, Po = — 2nj/(127r). However, in employing NNA, we promote Pq to the QCD value 
/3o = (33-2n^)/(127r). 

An alternative to the NNA procedure is to compute the vacuum polarization, including 
both the quark and gluon contributions, in the background-field gauge. The motivation 
for this approach is that, in the background-field gauge, the logarithmic dependence on the 
renormalization scale fi is contained entirely in the vacuum-polarization diagrams |T^. In 
the background-field gauge, the vacuum-polarization contribution still takes the form ( p. 161) . 
The coefficients a and b are computed in Appendix ^ and are given in Eq. 



B. Decomposition into Real and Virtual Contributions 



It is convenient to separate R^'^^ into real-real, real-virtual, and virtual- virtual contribu- 
tions. First, we make use of the identity 



lm{AB) = (lmA)B* + A(lmB) 

to identify the real and virtual parts of a single propagator: 

1 



Im[4m^ir('™)| 



X] 



X — le 



-Im [-in(x)]" - 7r5(x)Hn(0)]' 



Then, we substitute this result into Eq. ( p.l4|) to obtain 



(2.18) 



(2.19) 
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i?^"" = RoiG2R + 2GvG,R + Gl), (2.20) 

where the three terms in parentheses correspond, respectively, to the real-real, real-virtual, 
and virtual- virtual contributions. 

0.n - tJ^^C ^^n^,ym^)lt\vm -V^-V^) (2.21) 

is the contribution to the amplitude from two chains of real vacuum-polarization bubbles. 

/{j'") (x) = -2 Im [-mix]]"" (2.22) 

is the mth-order contribution to the amplitude for a single real bubble chain with virtuality 

X. 

GiR=tl'^fi^^O)I^R\^) (2-23) 

is the amplitude for a single real bubble chain, integrated over its virtuality and weighted 
with the heavy-quark and phase-space factor, with the virtuality of the other chain set to 
zero. 

oo 

Gv=T. (2-24) 

is the amplitude for a single chain of virtual vacuum-polarization insertions. In the last term 
of Eq. (|]2|), we have used the fact that /(O, 0) = 1. 

C. Manifestly Infrared-Finite Expressions 



In computing the expression for i?^"^ in Eq. (^M ), it is useful to re-arrange the individual 



contributions to make the expression manifestly infrared finite. First, we define the quantity 

Gi = GiR + Gv. (2.25) 

Since Gm and Gy are the real and virtual contributions to a single vacuum-polarization- 
bubble chain, the KLN theorem guarantees that Gi is infrared finite. We verify the finiteness 
of Gi by explicit calculation in Appendix |B[ It is also useful to define 
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2— — ^IR 



E 

m,n=l 



"1 dx 
/o 2nx 



dy 



(2.26) 



27ry' 

Since the contributions from the regions of integration near a; = and/or y = are equal 
in the first and second terms in Eq. ( p.26| ), G2 is infrared finite. Then, using Eqs. ( p. 201) , 
( p.23| ) and ( |2.26| ), we can write 



R"^^ = Ro{Gi + G2) 



(2.27) 



Each term in Eq. ( p.27| ) is infrared finite. 

The quantity Gi is computed in Appendix |B|. The resuh is 



^ 1 , 7ra^/3o 
Gi = — arctan ■ 



00 1 
+ V / dxgj 



1- as Pod 



X) 



where 



gn{x) = hm —Im[— 211(0;)]" = — Im asl3o{d — Inx + in) 



and 



d = hm Re (6 — a) /e. 



In naive non-Abehanization, 



(2.28) 



(2.29) 



(2.30) 



d = ln[f?/ (Ami)] - G , 



(2.31) 



while in the background-field gauge 



d 



1 



67 5 3,^n ^, 5 ^ 

Uf (P - 1) G 

12 18 4^^ ^3 



+ In 



4m3 



(2.32) 



By re-arranging Eq. (|2.26|) , we can obtain an expression for G2 in which every integral 
is separately infrared finite, and, hence, we can take e to zero: 



13 



A JO X Jo V 

m,n=l " " y 

+ -gUx) C %n{y)[f{x, y)-il- x)(l - y)] 9il~V^~^)]. 
Jo X Jo y 



(2.33) 



Eqs. ( |2.27| ), ( ^^.281 ), and ( |2.33| ) are the basis for the calculations in this paper. 



Taking into account the /i dependence of d [Eq. ( p. 30 )], one can verify straightforwardly 
that and a^(/i)G'2(/^) are independent if the running coupling constant a<j(/x) 

satisfies the one-loop renormalizat ion-group equation 

{fid/dfi)a,{fi) = 2/?o«s(/i). (2.34) 

III. THE NONPERTURBATIVE REGIONS 

So far, we have treated the effects of vacuum-polarization insertions as if they were 
entirely perturbative in nature. However, we expect perturbation theory to break down in 
the regions in which the gluon virtualities x and y are near zero, and, for on-shell gluons, in 
the regions in which energies Iq and ko are near zero. 

A. Single Bubble Chain 

We begin by investigating the region of small virtuality for a single bubble chain. That 
is, we consider the region of small x in the second term of Gi in Eq. ( p.28| ), namely, 

E / dxgn{x) = - E 

Near a; = 0, the Inx in Gib becomes unbounded, and the series fails to converge. This 
means that, in principle, the result that one obtains by carrying out the integration before 
the summation may be different from the result that one obtains with the opposite order 
of operations. In fact, as we shall see, interchange of these operations does produce a finite 
difference in the result. 
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Gib = / dxgn{x) = —y2 / dxlm as l3o{d — In X + in) . (3.1) 
"'0 TT ~, Jo L J 



Since 

A 

dxln"x~n!, (3.2) 

one might expect the terms in the series in Eq. to exhibit factorial growth, arising from 
the region of integration near x = 0. In fact, a careful calculation shows that the factorial 
growth cancels when one takes the imaginary part. Factorial growth in a bubble-chain series 
is associated with a renormalon singularity in the Borel transform of the series. Explicit 
calculation of the Borel transform of Gib shows that the renormalon singularity in the Borel 
plane vanishes when one takes the imaginary part 0. 

It is useful, in understanding both the dependence on the order of operations and the 
absence of factorial growth, to re-write the expression for Gu as a contour integral: 

1 °° /■ 

Gib = 7r-T. [asl3o[d - H-z)] (3.3) 

n=l •'^ 

The contour C runs from z = 1 — ie to z = Oto z = 1+ie. Term by term in the series, we can, 
without encountering any singularities, deform the contour C into the contour C, which is 
a circle of unit radius, centered at the origin, and traversed in the clockwise direction from 
z = 1 — ie to z = 1 + ie. Along the new contour G', the magnitude of z is never small. 
Hence, the logarithm in Eq. ( |3.3| ) is bounded, and it is clear that there is no factorial growth 
in the series. Furthermore, the perturbation sum is now uniformly convergent, and we can 
interchange the summation and integration without affecting the result. Hence, 

Gib = 7^ [ dz ^ . (3.4) 

21X1 Jc 1 — asPoyd — m[—z)\ 

Suppose, on the other hand, that we carry out the perturbation summation before de- 
forming the contour. Then we have 



TT Jo 



G\^ = — I dxlm 



1 — asf3o{d — Inx + iTT 
If we deform the contour to G', then we encounter the Landau pole at 

1 



If 1 

2^iJc^^ I- as(3o[d - ln(-^)] ' ^^'^"^ 



zq = — exp 



asPo 



(3.6) 
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and we must include its residue in the result: 



2m Ic 1 - asf3o[d - In(-z)] ^ a,/3o' ^^'^^ 
Then, using Eqs. ( p^) and ( p.7|) , we have 

Gib = G[f^ V- (3-8) 

asPo 

The contribution from the residue at the Landau pole [the second term in Eq. (|3l8|)] is 
just the difference in results between carrying out the integration before the summation 
and carrying out the summation before the integration. It displays an essential singularity 
in as and has the characteristic form of a power correction associated with a renormalon 
singularity. However, as we have seen, in this case, there is no renormalon singularity 
associated with the Borel transform of the series. 

The ambiguity in the result for Gib raises a question as to which order of operations, if 
either, is correct. If we carry out the integration before the summation, then the contour- 
deformation argument shows that there is no contribution from the region of small virtuality. 
In this case, asymptotic freedom suggests that the perturbation series is reliable. On the 
other hand, if we carry out the summation first, then there is no reason to trust the resulting 
expression: for fixed x sufficiently small, we are outside the radius of convergence of the 
geometric series and outside the regime of asymptotic freedom. Indeed, the Landau pole in 
the summed series is most likely an artifact of perturbation theory that would not persist in 
the true solution of QCD and, hence, the low-virtuality contribution from the residue at the 
pole is spurious.]^ Therefore, we conclude that the correct order of operations is to perform 
the integration first and then the summation. 



^Note that the situation is less clear in the present case than, for example, in the analysis of the 



r hadronic width, for which similar contour deformation arguments have been made jl3]. The 
r width possesses a spectral (Kallen-Lehmann) representation, which gives one some information 
about its analytic properties. The gluon propagator, being a colored object, possesses no spectral 
representation. 
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We note that the procedure of integrating over the gluon virtuahties before carrying out 
the perturbation summation is equivalent to the standard Borel-transform method for 
evaluating the contribution of a bubble chain. In that method, one takes the Borel transform 
of the series, sums the transformed series, carries out the integration over the phase space, 
and, finally, takes the inverse Borel transform. Since the initial Borel transform renders the 
series convergent for all virtuahties and absolutely integrable, one can interchange the order 
of operations, carrying out the integration first and then summing the Borel-transformed 
series. Furthermore, term by term in the series, the integration and Borel-transformation 
commute. Therefore, the standard Borel-transform method is equivalent to carrying out the 
integration, taking the Borel transform, summing the series, and then taking the inverse 
transform. This, by the definition of the inverse transform, is equivalent to carrying out the 
integration and then summing the perturbation series. 

In Appendix p, we have computed the first term of Gi in Eq. ( |2.28|) by performing the 
integration over x before the summation of the perturbation series. In this case, that is the 
only sensible order of operations: the perturbation series is not convergent for small values 
of the dimensional regulator e until one has canceled the infrared divergences by carrying out 
the integration over x in the real contribution and adding the real and virtual contributions 
order by order. 

We must also consider the possibility that, when x (y) is sufficiently close to unity, the 
energy /q (/^o) of the on-shell gluon may be too small to lie within the perturbative regime. 
(Recall that, although we have written Eq. ( |3.1| ), for compactness, only in terms of x, it 
derives from expressions that are symmetric in x and y.) Let us, for the moment, set aside 
the contributions from the nonperturbative regions of small /q (^o) by imposing a restriction 
on the region of integration: 

x,y < (l-VS)^. (3.9) 

Here, 5 is a positive parameter that we take to be much less than unity, but not so small 
that Am^S is outside the perturbative regime. Using Eq. ( p.9|) , we see that the restriction 
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9]) implies that 

ll> ml{2V5 - 5f ^ Amid, (3.10a) 
kl> ml{2VS - 6f ^ Amis. (3.10b) 



As we shall see in Sec. |III C| , the nonperturbative regions that we have set aside in imposing 
this restriction yield contributions of the form of one-loop corrections to the NRQCD matrix 
element of a color-octet 4-fermion operator. 

B. Two Bubble Chains 



The quantity G2 [Eq. ( |2.33| )] arises the from real-real contribution of the bubble chains. 



Hence, both virtualities x and y are nonzero. Furthermore, there are constraints on x and 
y. The first term in braces is subject to the constraint ^/x + y/y > 1; the second term in 
braces is subject to the constraint ^/x + y/y < 1. Owing to these constraints, one cannot 
immediately apply the contour-deformation argument of the last subsection to avoid the 
regions of small virtuality. For the first term in braces, when x is near unity, y ranges from 
unity almost to zero [(1 — -^/a;)^]- O'^s could write the integral over y as a. contour integral 
in the complex plane. However, the lower endpoints are tied down at ?/ = (1 — -^/a;)^ ± ie, 
and, hence, one cannot deform the contour out of the small-?/ region. For the second term 
in braces, when x is near unity, then y ranges from zero up to a small number [(1 — -^/a;)^]- 
Again, one could write the integral over y as a contour integral in the complex plane, but 
owing to the fixed endpoints at y = (1 — y^)^ i cannot deform the contour out of 

the small-y region. We deal with this difficulty by setting aside, temporarily, parts of the 
region of integration by imposing, again, the restriction (|3.9| ). 

For the first term in braces in Eq. ( p.33|) , the restriction (|3.9| ), combined with the con- 



straint y/x + y/y > 1, guarantees that 

x,y>S. (3.11) 
Hence, the x and y integrations lie entirely within the perturbative regime. 



For the second term in braces in Eq. ( |2.33|) , we may use the restriction ( p.9| ) to decompose 
the region of integration as follows: 

J dx J dy 9(1 — — y/y) = J dx J dy9{l~y'x — ^) 

+ dx dy + dx dy + dx dy . (3.12) 

Js Jo Jo Js Jo Jo 

In the first term on the right side of Eq. (p.l2|) , both integrations lie entirely within the 
perturbative regime. In the second term on the right side of Eq. ( p.l2| ), the x integration 
lies entirely within the perturbative regime, while the y integration does not. However, 
we can write the y integral, as in the preceding subsection, as an integral over a contour 
running from 5 — ieio zero to 5 + ie. We can then deform this contour into a circle of radius 
5 centered at the origin that is traversed in the clockwise direction from 5 — it to 5 + it. It is 
easy to see that one encounters no singularities in the integrand in performing this contour 
deformation. Hence, we conclude that there are no contributions to the second term on 
the right side of Eq. ( |3.12| ) from the region of small y. Similarly, for the third term on the 
right side of Eq. ( p. 121 ), we can write that the x integral as a contour integral and deform 
the contour out of the region of small x. Finally, for the fourth term on the right side of 
Eq. ( p.l2| ), we can write both the x and y integrals as contour integrals and deform both 
contours out of the region of small virtuality. Therefore, we conclude that, once we have 
imposed the restriction (|3.9|) , the contributions to the second term in braces in Eq. (p.33|) 
lie entirely within the perturbative regime. 

Again, we shall find in Sec. [HI C| that the contributions that we have set aside by virtue 
of the restriction ( |3.9| ) correspond to one-loop corrections to the NRQCD matrix element 
of a color-octet 4-fermion operator. We note that the contributions to G2 that we have 
set aside, when treated in perturbation theory, yield renormalon singularities. Such sin- 
gularities are a signal that the contribution cannot be computed reliably in perturbation 
theory. In Appendix |C|, we compute the leading renormalon singularity that arises from 
these contributions to 6*2- 
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C. Contribution of the Color-Octet NRQCD Matrix Element 



Now let us discuss the nonperturbative contributions that we have set aside by virtue of 
the restriction on the region of integration ( |3.9| ). For the contributions that have been set 
aside, we have either x > (1 — v^)^ or ?/ > (1 — v^)^, but not both, since this would violate 
the kinematic limit a/x + ^/y < 1. Hence, one of the final-state bubble chains is highly 
virtual compared to the momentum of the other bubble chain or the relative momentum 
of the incoming heavy QQ pair. Therefore, we may approximate these contributions by 
shrinking the highly virtual bubble chain to a point, replacing the bubble chain and its QQ 
interaction vertices with a four-fermion interaction. Corrections to this approximation may 
be taken into account by including four-fermion interactions containing derivatives with 
respect to the heavy-quark separations. It follows that the contributions that have been 
set aside are actually one-bubble- chain corrections to the matrix element of a four-fermion 
operator in the rjc state. 

In NRQCD at leading order in f , the four-fermion operator that arises when one shrinks 
a bubble chain to a point is the color-octet, spin-triplet operator 

Os{^Si) = ^^a'T^xX^f^'T^^lj . (3.13) 

We now verify that the one-bubble-chain corrections to the contributions of this operator to 
R do indeed reproduce, at leading order in f , the contributions that we have set aside. 

According to the NRQCD factorization formula [jl[ the contribution of this operator to 
the decay rate factors: 

(rye|08('^i)|^c), (3.14) 

where Im [/8(^S'i)/m^] is the short- distance coefficient. By matching Eq. ( |3.14| ) to the imag- 
inary part of the amplitude in full QCD for QQ QQ via a single bubble chain, one finds, 
in the bubble-chain approximation, that 

= -2gHui[K{l)], (3.15) 



2Im 



/8(^^l) 



2Im 



/8(^^l 



mt 
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where K{x) is defined in Eq. ( |2.7|) . 

The color-octet matrix element is fully determined by the NRQCD effective Lagrangian. 
It is nonperturbative in nature and must be determined phenomenologically or through 
methods such as lattice QCD. However, our aim here is to compute the perturbative contri- 
bution to this matrix element in the bubble-chain approximation so that we can compare it 
with the part of our perturbative calculation of i?^"^ that we have set aside. 

In NRQCD, the one-bubble-chain correction to the operator matrix element arises, in 
leading order in the heavy-quark velocity v, through the insertion of a ■ B operators into 
any two of the four quark lines. In the Feynman gauge, this contribution can be written as 

\ 2 



I 



-A;o-kV(2mc) +ie^ 



{H\0^{^S^)\H) 



^ OAT ^ \ ^0^^4 



'^Nc^J (27r)4 ^ ^ \^_A;o-kV(2m,) + zey 2ml 
x{H\0,CSo)\H). (3.16) 

The superscript A is a reminder that the NRQCD effective theory requires an ultraviolet 
cutoff. H denotes a hadronic state. 

The ko integral can be re-written by making use of the techniques of contour integra- 
tion. Closing the ko contour in the lower half plane, we pick up the residues at the poles 
and the discontinuities across the cuts in K^"^\x). The contribution of these residues and 
discontinuities can be written in terms of the imaginary part of K^'^\x). Then we have 



2 



'C'V m=0 



(27r)4 ' ' ^' y-ko-ky{2mc)+te^ 



x{H\0,CSo)\H). (3.17) 

Using Eq. (^.91), we can replace integration variables ko and k with x and y. We note 
that 

dkod\k\^ = imldxdy. (3.18) 
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Making this change of variables in Eq. ( p.l7| ) and multiplying by the short-distance quantity 
in Eq. ( |3.15| ), we obtain the one-bubble-chain correction to the rate: 

2TialCF ^ dx 



E 



Nc-ml ^0^0 27r 



2Im 



4mlK{l) 2lm 4mlK^"'\x) f{x,y){H\OiCSo)\H), 



(3.19) 



where f{x, y) is defined in Eq. ( p.l5| ), and we have used the approximate relation 



k2 , 

ko + - — ^ko- 



2mr 



2mr 



(3.20) 



which holds for ko/mc << 1. Dividing the rate in Eq. ( |3.19| ) by the rate into two photons 
in Eq. ( |2.13| ), we obtain the contribution of the color-octet matrix element to R^'^^: 



m=o^o 27r 27r 



4mlK{l)] 2lmkmlK^'^\y)] f{x,y). 



Now we choose the cutoff A on the matrix element to be the hard cutoff 



(3.21) 



x< S; 

y> (1 - vsy. 

From Eq. ( p.9|) , we see that this choice of cutoff corresponds to 

k^< [A;V(4me) + 2m^V6 - mj]'^ - fc^ 



(3.22) 



(3.23) 



Then, symmetrizing Eq. ( 3.21| ) in x and y, we find that it is precisely the contribution to 
Eq. ( ^.14| ) that is excluded by the restriction (^.9|). 

We conclude that the contributions that we have set aside by virtue of the restriction ( p.9| ) 
correspond to contributions to the decay rate that are proportional to the matrix element 
of the color-octet operator Os{^Si). In deriving Eq. (|3.21| ), we have dropped contributions 
of higher order in v. These could be taken into account by including the contributions of 
higher-dimensional color-octet operators. 
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We note that the use of a hard cutoff in the color-octet matrix element, as opposed to 
dimensional regularization, is essential in order to contain the small- virtuality contributions 
entirely within the matrix element 0. These contributions have an ultraviolet-divergent 
power behavior. In dimensional regularization, power-divergent contributions are set to zero, 
and, therefore, the small-virtuality contributions would be excluded from a dimensionally- 
regulated matrix element. 

Using Eqs. (|3.14| ) and ( p.l5| ), along with Eqs. (|L^ ) and ( p. 131) , we see that the contribu- 
tion of the color-octet operator Os{^Si) to R is given by 



Rs = -2 Im [AmlKil)]Ro ' X /io • (3-24) 









UsCf 




\ric) 



We note that, in the limit in which the dimensional regulator e is taken to zero. 



Im [AmlK{l)] = „, , r° , . . (3.25) 



In this paper, we will not compute Rs, but, instead, will treat it as an uncertainty in the 
calculation. The ratio of matrix elements in Eq. ( ^.241 ) may be estimated by making use of 
the velocity-scaling rules for NRQCD |I| . If we adopt the suggestion of Petrelli et al. []T1 



and include, in addition to the factors of the velocity f, a factor 1/{2N^ in the estimate of 
the color-singlet matrix element, then we obtain 

,3 



~ (3.26) 



2iV, 

A one-loop perturbative estimate gives a result that is numerically very close to that in 
Eq. (H): 

(3-27) 







{Vc\Oi 


?S,)\r,^) 



In making numerical estimates of the size of i?8, we will use the expression ( p. 27 



The value of the color-octet matrix element depends on the cutoff 5. In making the 
estimates ( pj.26D and ( |3.27| ), one assumes that the cutoff, in terms of momentum, is of 



the order of the dynamical scale mv. On the other hand, we would like to choose 5 so 
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that nonperturbative contributions arising from momenta less than v are contained in the 
color-octet operator matrix element, while contributions arising from larger momenta are 
contained in the short-distance coefficients. We will assume that a momentum cutoff of 
order mv is consistent with these desiderata. From Eq. ( p.23| ), we see that 



kVm^ < (2v^-5)2. (3.28) 

The choice S = 0.1 yields k^/m^ < 0.28. For charmonium, 0.28 is very close to the average 
value of w^. Therefore, we will use the value 6 = 0.1 in our calculations. 

IV. EFFICIENT METHOD FOR NUMERICAL COMPUTATION 

The perturbation series in Eq. ( p.28| ) is convergent. Once we have imposed the restriction 
( p. 91 ), the series in the first and second terms of Eq. ( |2.33| ) are also convergent. However, for 
reasonable choices of 6, the series converge very slowly. For example, for the second term 
in Eq. ( p.33| ) with 6 = 0.1 and = 0.247, it is necessary to compute through 12th order 



in as in order to achieve an accuracy of about 3%. Furthermore, at high orders in as, the 
integrand becomes strongly oscillatory, and the integral is difficult to compute numerically. 
It would be far more efficient if we could carry out the perturbation summation before 
integrating over x and y. 



We have already seen in Sec. [Ill A| that, in the case of Gi, we can carry out the perturba- 
tion sum first, provided that we compensate by adding a contribution that is proportional 
to the residue at the Landau pole. Using Eqs. (|3?T|) , (|375| ) and (|3.8|) , we can write 

dxgnix)= dxg{x) — , (4.1) 

n=i "'0 Jo asPo 



where 



9{x) = £ g^x) = , ""f" , , , . (4.2) 



n=l 



[1 - asf3o{d -\nx)] + (asPoir) 



Then, from Eq. ( |2.28| ), we have 
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G,(^) = ^arcta„-i^^4T;+/' <ix ^ + T^'' dxsW - (4.3) 

vrasPo I — as Pod J{i-^/5)^ x Jo asPo 



The second term in Eq. ( |4.^ ) arises from applying the cutoff ( ^.9| ) to GiRb [Eq- 

The analysis of Eq. ( p.33|) is somewhat more involved. For the first term in Eq. ( p. 331) 



the restriction on the region of integration ( p.9|) , together with the constraint ^/x + ^ < 1, 
guarantees that x,y > 6. Then the sums over n and m are absolutely convergent, provided 
that 6 is sufficiently large that 

\d -\nS + i7r\asPo < I. (4.4) 



For as = 0.247, Eq. (|4j) requires that 6 be greater than 0.0482 for the NNA value of d 
[Eq. ( p.31|) ] and greater than 0.0752 for the background-field-gauge value of d with ^ = 1 



[Eq. (|2.32| )]. Then, the convergence of the series and the absolute integrability of the sums 
allow us to interchange the summations and the integrations. 

For the second term in Eq. ( |2.33| ), we decompose the region of integration according 
to Eq. ( p.l2[ ). For the first term on the right side of Eq. ( |3.12| ), the perturbation series 



are absolutely convergent, provided that Eq. ( |4.4| ) is satisfied, and we can interchange the 
summations and the integrations. For the second term on the right side of Eq. ( |3.12| ), the x 
integration lies within the region of absolute convergence of the series, but the y integration 



does not. However, one can apply the contour-deformation analysis of Sec. [Ill A| to the y 
integral, deforming the contour to a circle of radius 6 centered at the origin, so that the 
entire contour lies within the region of absolute convergence. [The additional factors in the 
integrand of the second term in Eq. ( |2.33| ) contain no singularities within that circle.] It 
follows that we can interchange the y integration with the summations, provided that we add 
a compensating term that is proportional to the residue at the Landau pole. Similarly, for 
the third term on the right side of Eq. ( |3.12| ), we can interchange the integrations with the 
summations, provided that we add a compensating term to the x integral that is proportional 
to the residue at the Landau pole. For the fourth term on the right side of Eq. ( |3.12| ), we 
must add Landau-pole-residue terms to both the x and y integrals in order to interchange 
the integrations with the summations. Then, the various terms can be recombined to obtain 
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— (1 -x)g{x 

X 

(1-^)' dx , 
— 9\x 



dy 



y 



X 



-9{y) [f{x, y)-{l- - y)] e{l-^- v^) 
y 



2zo /-(i-v^)' dx 
^/ —9{x)[f[x,ZQ)-[l-x)[l-ZQ)\ 

ftsPo "'0 X 

+ ^Af{zo,z,)-{l-z,f]. 



(4.5) 



We have verified that the expressions (|4.3| ) and (|4.5| ) agree, within truncation errors, 
with numerical evaluations of the perturbation sums for Gi and G2 through 12th order in 
a,c. 



V. RESULTS 



In order to combine our result from the resummation of vacuum-polarization bubbles 
with the computation through order a^, we must subtract the part of the contribution 
through relative order «<, that is contained in the resummation. Expanding Gi [Eq. ( ^.28[ )] 
through order a^, we obtain 



Gi = 1 + a,/?o(l + d) + .... 



(5.1) 



The perturbation series for G2 [Eq. ( |2.33| )] begins at order a^. Therefore, using Eq. ( p. 271) , 
we see that we must subtract Rq[1 + 2as/?o(l + d)]. We conclude that the rate i?, including 
the exact next-to-leading-order contribution and bubble resummation of the higher-order 
contributions, is given by 



6 



TT 



+2a^(/i)/?o 



In 



4m? 



1 — d{jj) 



For the case of NNA resummation, we have 



(5.2) 



(5.3) 
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while for the case of background-field-gauge resummation, we have 

137r2 
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3 2^^ 



TT 



-2a^(/i)/5o- 



TT 



(5.4) 



where ^ is the gauge parameter for fields internal to loops. (.^ = 1 in the Feynman gauge.) 

We now evaluate R^^^, taking Uf = 3, choosing fi = 2mc, with rric the charm-quark pole 
mass, and using af\2mc) = 0.247 ± 0.012.0 We use Eqs. and (gj) to compute Gi and 
G2, respectively, with the cutoff 6 = 0.1. 

In the case of the NNA resummation, we obtain 



^NNA 



R 



NNA 



-0.64, 



(3.01 ±0.30 ±0.34) X 10^ 



(5.5) 



where the first uncertainty in R^^^ comes from the uncertainty in as{2mc), and the second 
uncertainty comes from the estimate of the uncalculated color-octet contribution [Eqs. ( p. 241) , 
O), and (13:271 )], with v = 0.3. 

In the case of background-field-gauge resummation, we take the gauge parameter ^ to 
be equal to unity, since this choice minimizes the size of the residual relative-order-a^ con- 
tribution. Then we obtain 



BFG 



1.88, 

r, - -1.04, 

^BFG ^ (^3_2Q ± 0.31 ± 0.47) X 10^ 



^BFG 



(5.6) 



^We have obtained this value by evolving with 3-loop accuracy from as{mz) = 0.118 ± 0.002, 
where mz = 91.188 GeV is the Z-boson mass, and using for the MS charm-quark and bottom- 
quark masses = 1.25 ±0.10 GeV and = 4.2 ±0.2 GeV, respectively The uncertainty 
in as(2mc) includes the uncertainty in as{Mz) and the uncertainty in the value of rric- 
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where the uncertainties in R^^^ are as in BJ^^^. 

There are additional uncertainties in our results that arise from the dependences on the 
renormalization scale fi. We expect the scale dependences in the resummed expressions to be 
considerably less than in the NLO expression. As we have already mentioned, -Ro(/^) [Gl{f^) + 
6*2 (/w)] is invariant with respect to changes of scale at the level of one-loop running of as{fi). 
The remaining, unresummed terms in R^'^'^ are a small fraction of the total expression and 
yield a correspondingly small dependence on the scale. If we evolve from as{2mc) = 0.247 
using the one-loop beta function, then we obtain as{mc) = 0.327. Taking this value of 
as{mc) and fi = rric, we obtain a 10.6% increase in the value of R^^^ and a 1.8% increase in 
the value of R^^^. At this value of ^^(mc), the inequality ([4.4|) is no longer satisfied for the 
background-field-gauge resummation, and the perturbation expansion no longer converges. 
For purposes of estimating the scale dependence, we assume that Gi and G2 are given by 
Eqs. ( |4.3| ) and ( [4.5| ) and, therefore, that i?o(/i)[Gf (/i) + G2{fi)] is scale invariant. 

Under three-loop evolution of ^^(yu), Ro{fi)[Gl{fi) + G'2(/^)] is no longer scale invariant, 
and, hence, the scale dependence is larger. Three-loop evolution from as{2mc) = 0.247 
yields as{mc) = 0.350. At this value of as{mc), the inequality ( ^^ ) is not satisfied in either 
the NNA or background-field-gauge resummations, and we again assume, for purposes of 
estimating the scale dependence, that Gi and G2 are given by Eqs. (|47^ ) and ([4.5|) . Taking 
the three- loop-evolved value of as{mc) and fi = rric, we obtain a 24.5% increase in the 
value of R^^^ and an 11.9% increase in the value of R^^'^. One could hope to reduce 
these uncertainties further through a resummation of two-loop bubble contributions. Still, 
they are already a significant improvement over the factor of 2.4 uncertainty from scale 
dependence in the NLO result. 

VI. SUMMARY AND DISCUSSION 

We have studied the ratio R of the hadronic and electronic decay rates for the heavy 
quarkonium rjc- In R, the dependences of the rates on NRQCD matrix elements cancel at 
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leading order in v. The relative-order f ^ corrections also cancel. Consequently, renormalon 
ambiguities associated with initial-state virtual gluons cancel at this level of accuracy in 

The perturbative corrections to R that are associated with the running of as account for 
most of the large one-loop correction to R. We have carried out a resummation of the leading 
corrections of this type, which arise from vacuum-polarization insertions in the final-state 
gluon legs. Our results are in good agreement with the experimental value of R. 

There is some arbitrariness in the choice of the resummed quantity. In this paper, 
we have used the NNA approach 0, in which one resums, in nth order, all contributions 
proportional to (as/^o)". In implementing this method, one resums chains of quark-loop 
vacuum-polarization bubbles in gluon legs and then accounts for gluon contributions by 
promoting the quark contribution to (3o to the full QCD /?o- An alternative to the NNA 
approach is to resum chains of quark-loop and gluon- loop vacuum-polarization bubbles in the 



background-field gauge |T0[. This method is motivated by the fact that, in the background- 



field gauge, all of the running of is contained in the vacuum-polarization diagrams [IT2 



The two approaches give numerical results that agree at the level of the other uncertainties 
in the calculations. 

The sum of the vacuum-polarization-bubble terms is independent of the renormalizat ion- 
scale, at the level of the one- loop running of a^. Hence, the renormalization-scale dependence 
of our result is much milder than that of the unresummed one- loop result. However, because 
of the arbitrariness of the resummation procedure, it would not be correct to conclude that 
the reduction in the theoretical uncertainty is proportional to the reduction in the scale 
dependence. 

The final-state phase space includes integrations over the virtualities m'^x and m'^y of 
the final-state gluons. Near zero virtuality, one is in a region in which perturbation theory 
breaks down and nonperturbative contributions may be important. We have identified this 
region of low virtuality with contributions to the NRQCD matrix element of a color-octet, 
^Si four-fermion operator. The contribution of this matrix element is of relative order v^, 
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and it is treated as an uncertainty in our calculation. 

A key ingredient in our analysis is the demonstration that there are no contributions 
in which both gluons have small virtualities or in which one gluon is exactly on shell and 
the other has a small virtuality. Such contributions, if present, would not correspond to an 
NRQCD operator matrix element and would be inconsistent with the NRQCD factorization 
program. We demonstrate the absence of such contributions by writing the integrals over 
the virtuahties of the off-shell gluons as contour integrals and deforming the contours out of 
the regions of small virtuality. 

One can understand this contour-deformation argument in terms of the uncertainty prin- 
ciple. In the bubble-chain picture, away from the kinematic limits of maximum gluon vir- 
tuality, the quarkonium decay produces two energetic gluon jets, each with energy of order 
rric- Consider the case in which both gluons are virtual. That is, each gluon jet contains two 
subjets that correspond to the constituents of a vacuum-polarization bubble. The mass of 
each jet ranges up to a value of order rric. Now suppose that the jets are observed inclusively, 
that is, without regard to their masses. Provided that neither jet has a mass near the kine- 
matic limit, the virtual-gluon masses are independently uncertain by amounts of order rric 
and, according to the uncertainty principle, the virtual gluons propagate over distances of 
order l/rric- Hence, the NRQCD four-fermion annihilation vertex is well localized, in accor- 
dance with the NRQCD factorization formula. On the other hand, if one restricts the gluon 
virtualities {i.e., x and y) to be small, then the gluons propagate over long distances, of 
order l/{mc\/x) and l/{mc^/y). Such propagation would delocalize the annihilation vertex. 
Hence, contributions from small x and y must vanish in the inclusive rate. The absence of 
singularities in the complex virtuality plane near the vacuum-polarization cut tells us that 
there are no competing, nearly classical processes that would upset this uncertainty-principle 
analysis of the bubble chains in isolation. 

Note that this uncertainty-principle argument does not rule out contributions in which a 
gluon is exactly on the mass shell (unless the other gluon is virtual and at small virtuality). 
One should regard the on-shell gluon not as propagating an infinite distance, but, rather, 
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as following a nearly classical trajectory. Since the momentum and energy are uncertain by 
amounts of order rric, the positions and times in the trajectory are smeared by amounts of 
order l/rric, and one again reaches the conclusion that the NRQCD four-fermion annihilation 
vertex is localized to within l/rric . 

On the other hand, if the mass of one of the jets is near the kinematic limit, then the mass 
(and uncertainty in mass) of the other jet is constrained kinematically to be small, and the 
annihilation vertex becomes delocalized. As we have seen, this kinematic region, in which 
one virtuality is near the kinematic limit, yields long-distance (delocalized) contributions 
that are absorbed into color-octet operator matrix elements. 

An interesting, and somewhat unexpected, feature of the bubble-chain contributions is 
that their values depend on whether one carries out the integration over the gluon virtualities 
before or after carrying out the perturbation summation. The two orders of operations lead 
to results that differ by an amount that has the characteristic structure exp[— l/(as/3o)] of a 
renormalon ambiguity. However, this difference does not arise from a renormalon ambiguity 
since, in our calculations, all of the renormalon contributions are absorbed into the color- 
octet, ^Si matrix element. Rather, the difference arises from the failure of the perturbation 
series to converge in the region of small virtuality. The difference is associated with the 
presence of the perturbative Landau pole in that region. We have argued, again using 
contour integration, that one should integrate over the virtualities before carrying out the 
perturbation summation. Then one can deform contours out of the region of small virtuality 
and the resulting perturbation series is reliable. On the other hand, there is no justification 
for summing the perturbation series first: at sufficiently small virtuality, one is outside 
the radius of convergence, and the series is unreliable. We have argued that the order of 
operations in which one carries out the integration over the gluon virtualities before carrying 
out the perturbation summation is equivalent to the standard Borel-transform method for 
evaluating the contribution of a bubble chain. 

For as evaluated at the scale /x = 2mc and for typical values of the cutoff 6, the perturba- 
tion series, after integration over the gluon virtualities, converges very slowly. Furthermore, 
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the integrands in high-order terms oscillate strongly and are difficult to evaluate numerically. 
Therefore, we have devised a more efficient method for numerical computation in which one 
carries out the perturbation summation before the integrations over the gluon virtualities 
and then makes corrections to account for the interchange of the order of limits. 

Our numerical result for R is much smaller than that which one obtains by applying 
the BLM scale-setting method to the next-to-leading order calculation.^ This may seem 
surprising, since the BLM method is designed to recover the effects associated with the 
one- loop running of a^, to the extent that the can be taken into account through a change 
of renormalization scale. The effects of a change of scale in as can be expanded in a 
perturbation series. At one-loop accuracy in the running of a^, one obtains the geometric 



series 



«.(/^') = . , "^y, , ; . = f:[2a.(/.)/3oln(/.//x')]^ (6.1) 

It is illuminating to compare this series with the perturbation expansion for the contribution 
Gi{^) [Eq. ( |2.28| )] of a single bubble chain, integrated over the gluon virtuality. In the NNA 



case we have 



Gi{2m^) = 1 + 1.91a,(2mc) + 2.47a^(2mJ + 0.97a^(2mJ - 4.49a^(2mJ 

-11.76a^(2mc) + ■ ■ ■ . (6.2) 



We see that the individual terms in the expansion ( |6.2| ) are not well approximated by those 
of a geometric series. Therefore, we would not expect the BLM method to reproduce this 
result accurately. Of course, before integration over the gluon virtuality, the bubble chain 
is a geometric series [Eq. (|2.7| )], with a ratio between terms of —ill{x). However, in order 



one regards the BLM procedure as an approximate resummation of bubble contributions at 
one-loop accuracy in the bubble, then it would be appropriate to evolve Og from the nominal scale 
2mc to the BLM scale by using the one- loop (3 function. This leads to ii^^°(/iBLM) = 6.3 x 10^, 
which is smaller than the value, obtained by 3- loop evolution, in Eq. (|l.7| ). 
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for the bubble-chain contributions, integrated over phase space, to be close to a geometric 
series, we must have 

([n(x)n^(n(a:)r, (6.3) 

where the angular brackets denote averaging over phase space, with the other factors in 
the integrand for i? as a weight. The relation (|6.3| ) holds only if the integration over the 
virtuality x is highly peaked. Hence, we do not expect, in general, to obtain a geometric 
perturbation series. We note that there are generalizations of the BLM method that take 



into account higher-order deviations of the perturbation series from a geometric series |T5 
and that may reproduce the effects computed in this paper. 

Finally, we remark that the methods that we have developed in this paper can be applied 
to a variety of processes to resum large final-state corrections associated with the running of 
as- The methods are most directly applicable to computations in NRQCD, but with small 
modifications, they are also applicable to computations in full QCD. The principal change 
in the method that is required in the latter case is to absorb nonperturbative contributions 
into QCD matrix elements, rather than into NRQCD matrix elements. In general, there 
are contributions that are associated with the running of in initial-state virtual gluons, 
as well. Such initial-state contributions also contain nonperturbative pieces that must be 
absorbed into operator matrix elements. However, the bubble series begins at relative order 
a^, instead of relative order a^. 
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APPENDIX A: THE GLUON SELF-ENERGY IN THE BACKGROUND-FIELD 

GAUGE 



At one-loop level, a gluon self-energy receives contributions from gluon, ghost and quark 
bubbles. The value of a gluon bubble depends on the choice of gauge. In this Appendix, we 
describe the calculation of the gluon self-energy in the background-field gauge . 

In the background-field-gauge method, the gauge field is separated into an external field 
or background field, which appears only in external lines of a Feynman diagram, and an 
internal field, which appears only in loops of a Feynman diagram. The gauge of the internal 
field is fixed, but that of the external field is not. Hence, the action is explicitly gauge co- 
variant with respect to the external field, and the background-field-gauge method maintains 
explicit gauge covariance with respect to the external field in each step of a calculation. 
In the background-field gauge, all of the logarithmic dependence on the renormalization 
scale that is associated with the running of the coupling constant is contained in the gluon 



self-energy |12 



In the class of Lorentz-covariant background-field gauges, the gauge of the internal 
fields is specified by the parameter ^, with ^ = 1 corresponding to the Feynman gauge 



r^ . In the background-field gauges, the contribution of a gluon bubble to the gluon 



self-energy is given by 

Yfabgluonfi, c\ — ^ ^2 , ,2e racd rbdc f ^^V ^tiyiP^^iC) (\-[\ 

IV (^,4)--2^?s/i / / ] ^2^)^ fi^p + kf ' ^^^^ 

where 



X I (2p + k)^g^s + '2k^9uS - '^ksg^^ + [(^ + P)i9u& + Psgu-f] \ 



{p + kf 

A straightforward calculation, using Feynman parameter integrals, gives 



(A2) 
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X r(e)r(i-e) 



r(3-e) 2-^(1-6) 



r(4-2e) r(2-2e) 



-xe-i) 



(A3) 



The contribution of a ghost bubble to the gluon self-energy in the background- field gauge 



IS 



Y[abghost(i.\ _ „2 , 2e e acd ebdc f ^^P (2p + fc)^ (2p + fc) 



27T 



r(4 - 2e) 



(A4) 



The contribution of Uf massless quarks to the gluon self-energy is 

ab,uark(.. _ 2 2e^ (^.^l.. f d°P Tr[(^+ K)7m M 

{k}- g.njii irU ^ j p2 (p + k)2 

^ ^ TT \ — — le 



(A5) 



r(4 - 2e) ■ 

Adding (|A3|) , ( |A4|) and (|A5|) and subtracting the pole in e that corresponds to an ultra- 
violet divergence, along with an associated constant that depends on the renormalization 
scheme, we obtain the expression for the renormalized one-loop gluon self-energy: 



(A6) 



where 



r(i + e) 



r(2 - 2e) " ll" + r(4 - 2e) ~ 4^^^ ~^ 
67 5 3 



1 + 



—7 + In 47r + 



-mf — ze 



/?ovr V12 18 ' 4 



4m^ 



e + 0(e2) e'-, (A7a) 



a = l+ (-7 + ln47r + ^ + C)e + 0{e^) , 



(A7b) 



/3o = (33 — 2nj)/(127r), x = p^/(4m^), and C is a renormalization-scheme-dependent con- 
stant. In the MS scheme, C = —5/3. Now, in Eqs. (^) and ([A7|), e plays the role of an 
infrared regulator (e < 0). 
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APPENDIX B: COMPUTATION OF A SINGLE BUBBLE CHAIN 



In this Appendix, we compute the real and virtual contributions from a single bubble 
chain. We consider first the real contribution. From Eq. ( |2.23| ), we have 



G 



IR 



-Us (bx~ 



(Bl) 



where e < is an infrared regulator. The two terms in parentheses have different infrared 
behaviors. The contribution of the second term contains infrared divergences, which appear 
as poles in e. The contribution of the first term is free of infrared divergences, and, so, we 
can take e = in the integrand. Thus, the contribution of the first term can be written as 



G 



IRa 



— ^ / dxlm 



asPo {d — \iax + in) 



(B2) 



where d is defined in Eq. (|2.3CI| ). 

Now let's evaluate the second term in parentheses in Eq. 



T|). It can be written as 



G 



IRb 



vr ±r; Jo X 



n=l 

oo 



V^^Im 
If i^Im 



i+1 



[hy - a 

n— 1 



[-a) 

^by _ ay 



(B3) 



m=0 



where we have made a change of integration variable y = x ^, and, in the second line, we 
have used the fact that a is real. Carrying out the integration over y, we obtain 



G 



IRb 



m=l 



n=l 



1 g i-asM- 



TT 



n=l 

oo 



cn+1 



1 



a — b 



n+l 



a n + l 



n 

+ 0{e), 



(B4) 



n=l - 7ia,(3o^2^ \ - / 

where, in the last line, we have used the fact that the phase of b, exact to all orders in e, 
is ire. In the last line of Eq. ([B4D, the first term contains infrared divergences, and, since 
b — a ^ e, the second term is infrared finite. 
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From Eq. ( |2.24|) , we see that the virtual correction from a single bubble chain is 



G 



V 



E 

n=l 



(B5) 



The terms in Eq. ( [B5| ) with n > 1 exactly cancel the the infrared divergences in the first 
term of the last line of Eq. ( P^ ) . The n = term of Eq. ( |B5| ) and the remainder of Eq. 
combine to give 



GiRb + Gv 



-1 



-Im In 



1 - asPo 



b — a 



1 TTOsPo 

arctan ■ 



1 - asPod 



Combining Eqs. ( [B3| ) and (BU), we obtain 



Gl— GiRa 

1 



GiRb + Gy 



arctan 



1 — asPod vr 



H — E ' dxlm 



n=l 



as(3o{d — Inx + in) 



(B6) 



(B7) 



APPENDIX C: RENORMALON IN THE SHORT-DISTANCE COEFFICIENTS 

In this Appendix, we compute the leading renormalon singularities in the Borel trans- 
forms of the short- distance coefficients of We consider the expressions for Gi and G2, 
Eqs. ( |2.28|) and (|2.33| ), without the cutoff (|3.9|) . Once the cutoff ( p.9|) has been imposed, the 
renormalon singularities are absent in Gi and G2- 

The Borel transform of a function /(a^) that has a perturbation expansion 



(Cl) 



n=0 



is defined by 



B[f]{u) = ao5iu/Po) + J2 



(C2) 



ntl(n-l)! 

The renormalon singularities in the Borel transform are those that arise from the bubble- 
chain sum, and the leading singularity is the one that appears at the smallest value of 



u. 



^For reviews of the properties of renormalons and Borel transforms, see Refs. pi 
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Renormalon singularities can arise from the ultraviolet or infrared regions of momentum 
integrals. Since, in R, the virtualities x and y are bounded by unity, there are no ultraviolet 
renormalons. We have argued (in Sec. [Ill A| ) that Gi receives no contributions from the 
region of small virtuality. Hence, it contains no renormalons. We have also argued (in 
Sec. [Ill BD that G2 receives no contributions in which both x and y are small. Therefore, the 



renormalons in G2, can arise only from the regions x ^ and y ~ 1 or ?/ — and x ~ 1. 
We consider the region x — and y ~ 1 and multiply our result by two to account 



for the contribution of the second region. In this region, in the expression (|2.33|) , we can 
approximate gn{y) by Qni^)- Corrections to this yield behavior in x that is subleading as 
X 0. Furthermore, the 5'„(1) yield a convergent series, and, so, we may carry out the 
Borel transform term by term. From the relation 

B[asf]iu) = ^ r du' B[f]{u') , (C3) 

we see that the leading singularity in the Borel transform of G2 must come from the leading 
power of as in the terms Qni^)'- the integration on the right side of Eq. ( P3| ) weakens the 
singularity as the power of from the Qni^) increases. Therefore, we retain only the leading 
term gi{l) = a^/^o, to obtain 

Ib[G2]{u) -Po C-{l-x)B[asg]{x,u) [' ^-{l-y) 
2 Jo X J{i-y^Y y 

dx . ^ . . /-(i-v^)' dy 



+Po -B[asg]{x,u) ^ [/(^, ^) _ (1 _ ,:)(l _ ^)] , (C4) 

Jo X Jo y 



where 



B[asg]{x,u) = -Im (a;-«e'^"+^™ 



= x-^e""" sin(7rM) . (C5) 

An analysis of the integrations over y yields the leading contributions for x small. Inserting 
these results into Eq. (P^), we have 

1 dx dx 

-B[G2]{u) -Po —{l-x)B[asg]{x,u){2/3)x^/^ + po —B[asg]{x,u)3xlnx, (C6) 

^ J X J X 
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where the first and second terms on the right side of Eq. (|Cq ) correspond to the first and 
second terms, respectively, on the right side of Eq. ( |C4| ) . Then, a straightforward calculation 
gives the leading renormalon singularity for each of the two terms in Eq. 



B[G2]{u) ~ ^e'^" 

TT 



8 sin(7rn) sin(7ru) 
6- 



(C7) 



3 3-2M "(1-m)2_ 

where the first and second terms on the right side of Eq. (|C7| ) correspond to the first and 
second terms, respectively, on the right sides of Eqs. (|C4D and (|2.33| ). 
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